library(tidyr)
library(dplyr)
library(ggplot2)
library(ggpubr)

if (!file.exists("output")) {
  dir.create("output")
}

source("data_handling.R")

SBS_earlylate <- read.csv("ExtendedDataFigure_4_data.csv")

# calculate frequency of positive cases for either early or late signatures
signatures <- names(SBS_earlylate[-c(1, 2)])

frequency_sig <- NULL
for (sig in signatures) {
  n_pos <- SBS_earlylate %>%
    select("donor_id", "time", sig) %>%
    spread(time, sig) %>%
    filter(!(early == 0 & late == 0)) %>%
    nrow()

  freq <- data.frame(signature = sig, n_pos = n_pos, total = nrow(SBS_earlylate) / 2) %>%
    mutate(
      freq = round(n_pos / total * 100, 2),
      freq = paste0(n_pos, "/", total, " (", freq, ")"), .keep = "unused"
    )
  frequency_sig <- rbind(frequency_sig, freq)
}

# calculate mean and SD
mean_earlylate <- SBS_earlylate %>%
  group_by(time) %>%
  summarise_if(is.numeric, mean)
sd_earlylate <- SBS_earlylate %>%
  group_by(time) %>%
  summarise_if(is.numeric, sd)

distribution_sigs <- data.frame(
  signature = signatures,
  mean_early = as.numeric(mean_earlylate[1, -1]),
  sd_early = as.numeric(sd_earlylate[1, -1]),
  mean_late = as.numeric(mean_earlylate[2, -1]),
  sd_late = as.numeric(sd_earlylate[2, -1])
)

# Wilcoxon test
output_stats <- NULL
for (sig in signatures) {
  wilcox <- wilcox.test(get(sig) ~ time, data = SBS_earlylate, paired = T)
  stats <- data.frame(signature = sig, p_val = wilcox$p.value)
  output_stats <- rbind(output_stats, stats)
}

output_stats <- output_stats %>% mutate(p_adj = p.adjust(p_val, method = "BH"))

merge_table <- frequency_sig %>%
  left_join(distribution_sigs, by = "signature") %>%
  left_join(output_stats, by = "signature")

write.csv(merge_table, "output/Supplementary_Table_12.csv", row.names = F)
